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Abstract 

A fully adaptive methodology is developed for reducing the complexity of large dissipative sys- 
tems. This represents a significant step towards extracting essential physical knowledge from 
complex systems, by addressing the challenging problem of a minimal number of variables needed 
to exactly capture the system dynamics. Accurate reduced description is achieved, by construction 
of a hierarchy of slow invariant manifolds, with an embarrassingly simple implementation in any 
dimension. The method is validated with the auto-ignition of the hydrogen-air mixture where a 
reduction to a cascade of slow invariant manifolds is observed. 

PACS numbers: 47.11.-j, 05.20.Dd, 05.70.Ln 



*Elcctronic address: eliodoro.chiavazzo@polito.it 
^Electronic address: karlin@lav.ma vTethz.chl 



1 



I. INTRODUCTION 



Detailed reaction mechanisms typically serve as accurate models of dissipative complex 
systems with many interacting components: Biochemical processes in living cells and com- 
bustion phenomena are prototypical examples of such systems [l-^. Modern research has 
to cope with an increasing complexity mainly in two aspects: First, the number of degrees 
of freedom (scaling with the number of components) is tremendously large; second, complex 
system dynamics is characterized by a wide range of time-scales. For example, the usage of 
detailed reaction mechanisms in the reactive flow simulation soon becomes intractable even 
for supercomputers, particularly in the turbulent combustion of even "simplest" fuels such 
as hydrogen |4H6|. As a result, there is a strong demand for methodologies capable of both 
drastically reducing the description of complex systems with a large number of variables, 
and concurrently allowing physical insights to be gained. Modern automated approaches to 
model reduction are based on the notion of low dimensional manifold of the slow motions 
(slow invariant manifold - SIM - for short) in the phase-space describing the asymptotic 
system behavior. Although several methodologies have been suggested in the literature {zj], 
the construction of accurate reduced description remains a rather challenging task. In par- 
ticular, the evaluation of numerical SIM approximations in the phase-space is hindered by 
several difficulties as far as the choice of the manifold dimension is concerned, since the 
latter information is typically not known a priori. In addition, accurate simplification of 
complex multiscale systems often requires the construction of heterogeneous (variable di- 
mension) manifolds with the dimension d ranging from unity up to tens in different regions 
of the phase-space. To the best of our knowledge, at the present, fully adaptive model 
reduction methodologies capable to cope with the above issues are still missing. This re- 
search area is pretty active and much effort has been devoted to devising techniques with the 
above features. The intrinsic low dimensional manifold (ILDM) approach Q], the computa- 
tional singular perturbation (CSP) method and the minimal entropy production trajectory 
(MEPT) method 10J are only some representative examples. In addition, the minimal num- 
ber of reduced degrees of freedom underling the asymptotic dynamics of complex multiscale 
systems is still a debated issue ll|. In this respect, we notice that, though here we mainly 
focus on chemical kinetics, our results have direct implications on the study of the homoge- 
neous isotropic Boltzmann equation which has been stated a fundamental problem of Physics 
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12[. The latter investigation is beyond the scope of this work, however future works shall 



move in this directions, where we can take advantage of recently introduced models such as 



the one proposed in [13]. 



In the present work, we introduce a methodology which enables to cope with the accurate 
reduced description of large dissipative systems, where no a priori assumptions on the least 
number of fundamental (slow) variables are made. Toward this end, both global and local 
construction of slow invariant manifolds, with an embarrassingly simple implementation up 
to any dimension, is worked out. 

This paper is organized in sections as follows. In the section [Til we briefly review the 
governing equations for chemical kinetics. The problem of model reduction, as understood 
by the Method of Invariant Manifold (MIM), is discussed in the section HTT1 The Relaxation 
Redistribution Method (RRM) is introduced in the section IIV[ where both a global (section 
IIV A I) and a local (section IIV Aj) formulation are presented. The latter methodology is 



validated for a detailed chemical kinetics describing a reacting mixture of hydrogen and air 
in the section [V] Finally, conclusions are drawn in the section IVII 

II. DISSIPATIVE REACTION KINETICS 

In the present study, we assume that a complex dissipative dynamics is governed by an 
autonomous system in terms of the state fona phase space U with a unique steady state, 

Important example of (CQ) to be addressed below is the reaction kinetics where ip = 
(ipi, . . . ,ip n ) is a n-dimensional vector of concentrations of various species, while the vector 
field / is constructed according to a detailed reaction mechanism as described below. More 
specifically, in a closed reactive system, the complex reaction of n chemical species A±, A n 
and d elements can be represented by a (typically) large number r of elementary steps: 

n n 

J2 a siAi^J2^ si A h s = l,...,r, (2) 

i=l i=l 

where a S i and (3 si are the stoichiometric coefficients. The latter coefficients enable to define 
the three stoichiometric vectors: a s = (a s ±, a sn ), (3 S = (j3 s i, (3 sn ) and 7 S = (3 S — a s , 
where the index s runs over the r elementary reactions (T5]). For clarity, in the detailed 



reaction mechanism for air and hydrogen to be considered below [3[, s identifies any of the 
21 reactions in Tabled while the corresponding stoichiometric coefficients a S i and /3 S i indicate 
the number of molecules of species % in the reactants and products of reaction s, respectively. 
Production (or depletion) rates of chemical species can be conveniently expressed in terms 
of the differences: = (3 si — a si . 

Expressing the state in terms of the molar concentrations ip = (ci, ...,c n ) (ratios of the 
number of moles by the volume) , all chemical species evolve in time according to the mech- 
anism 02]): 

^ = X>^(^)' ( 3 ) 

where W s (ip) is the reaction rate function of the reaction s, which (usually) takes a polyno- 
mial form according to the mass action law: 

n n 

W s (VO = W+ (V, 9) - W- (V, 9) = k+ (9) J] cf - kj (9) J] cf , (4) 

i=l i=l 

with the reaction constants and k~ depending on the system temperature 9 according to 
the Arrhenius equation: 

k s (9) = A s 9 n *e- Ea * /ne , (5) 

where the quantities A s , n s , E as are fixed (and tabulated, see e.g. Table [Tj) and referred 
to as pre-exponential factor, temperature exponent, activation energy of the reaction s, 
respectively, while 1Z is the universal gas constant. Due to the principle of detailed balance, 
a relationship between the latter reaction constants (kf, k~) is established for each step s 
at the steady state: W£ = W~ . In general, the system ([3]) is to be solved in combination 
with an additional equation ruling the temperature evolution (energy equation). 

The concentration of the z-th chemical species can be also expressed in terms of the mass 
fraction = UiCi/p, while, in an adiabatic closed system, the temperature is computed by 
conserving the mixture-averaged enthalpy, which for ideal gases reads 

n 

h = ^2Yihi{e), (6) 

i=l 

where p, Ui and hi are the mixture density, the molecular weight and specific enthalpy 
(per unit mass) of species i, respectively. For the sake of completeness, we report here the 
closed dynamical system governing closed reactive ideal mixtures under fixed enthalpy h 
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and pressure p to be addressed below in section [V] 

[ dip/dt = E lsW s {iP, 6) = (dd/dt, dc n /dt) 

I 8=1 n (7) 

|^ de /dt = -±f J h i (9)Y t 

where C p denotes the mixture-averaged specific heat under fixed pressure, while specific 
enthalpy hi (6) for any species i can be computed using fflUj) . Molar concentrations q are 
linked to mass fractions Yj as q = p(Yi/ui) / (lZ9 E™ Yj/ujjj , while the mass fraction rate 
Yi reads as follows: Yi = Uif)~ 1 dci/dt, where Ui is the molecular weight of species i. We 
notice that the second equation in ([7]) stipulates the conservation of h, thus it represent an 
alternative way of imposing Constance of (JH])- 

Finally, due to the conservation of elements, in a closed reactor, d linear combinations of 
the species concentrations (expressing the number of moles of each element) remain constant 
during the system evolution in time: 

Cip = const, (8) 

where C is a d x n fixed matrix. 

Remark-Having in mind dissipative multiscale dynamics such as chemical and physical 
kinetics, here we focus on systems ([T|) with a single steady state. Hence, the Relaxation 
Redistribution Method (RRM) introduced below in section IIVI has been tested for those 
cases so far. We stress however that, for deriving the RRM approach, no assumptions are 
made concerning the number of steady state points of (CQ). Thus, implementations of the 
RRM to different dynamics shall be presented in future publications. 

A. Thermodynamic Lyapunov function 

Due to the second law of thermodynamics, the kinetic equations ([3]) are equipped with 
a global thermodynamic Lyapunov function G(ip). In other words, the time derivative of 
the above state function is non-positive in the whole phase-space, G (ip) < 0, with the 
equality holding at steady state. For instance, in an adiabatic reactor with fixed pressure 
p and enthalpy h, the specific mixture-averaged entropy s (in mass units) monotonically 
increases in time starting from any non-equilibrium initial condition: hence the function 
G = — s decreases during the dynamics. For ideal gas mixtures, a Lyapunov function G of 
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the system ([3]) takes the explicit form: 

n 

G = -s = - x i & W - KlnXi - Kin (p/pref)} /W, (9) 

i=l 

where X{ — c^j Ci anc ^ Si denote the mole fraction and the specific entropy of species 
i, respectively, 1Z is the universal gas constant, p re f a reference pressure and W the mean 
molecular weight. For numerical purposes, the properties of the i-th species, hi and Sj, can 
be expressed in terms of the temperature, 9, and a set of tabulated coefficients as follows 

hi (9) = K9 (an + ^9 + ^f9 2 + ^f9 3 + 2f # 4 + *f) , 
Si (9) = K (a u \n9 + a i2 9 + <f9 2 + <f9 3 + 9 A + a i7 ) . 

III. THE FILM EQUATION OF DYNAMICS 

If the number of degrees of freedom n is large, one may seek a reduced description with a 
smaller number of variables q <n. A consistent approach to model reduction is provided by 
the Method of Invariant Manifold (MIM) whose brief review is in order. Interested reader 



is delegated to the work [15] for further details. 

In MIM, the problem of model reduction is identified with the construction of a slow 
invariant manifold (SIM) fisiM, whose dimension q is the number of the essential (macro- 
scopic) variables which parameterize the SIM. As sketched in the cartoon in Fig. Wfa), the 
above method is based on the idea that the macroscopic slow dynamics of a complex system 
occurs along the SIM (invariance), once an initial fast relaxation toward the SIM has taken 
place. Let a manifold Q (not necessarily a SIM) be embedded in the phase space U and 
defined by a function Q = which maps a macroscopic variables space 5 into U. Intro- 
ducing a projector P onto the tangent space T of a manifold Q, the reduced dynamics on 
it is defined by the projection Pf(Q) e T (see Fig. [lb)). A manifold Q is termed invariant 
(but not necessarily slow) if the vector field / is tangent to the manifold at every point: 

/MO) - Pfm)) = o, £ e s. 

While the notion of a manifold's invariance is relatively straightforward, a definition of 
slowness is more delicate as it necessarily compares a (faster) approach towards the SIM 
with a (slower) motion along SIM. In MIM, slowness is understood as stability, and SIM is a 
stable stationary solution ^sim(0 °f ^ ne following film equation of dynamics defined on the 
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FIG. 1: (Color online), a) Model reduction techniques assume the following idea: After a fast 
initial transient at the time instants t < to, the (slow) dynamics of a complex system takes place 
along a slow invariant manifold (SIM) on the phase-space U at any future time t > to (invariance) 
toward the steady state, b) The definition of a projector P onto the tangent space T introduces a 
decomposition of slow and fast motions of the field /. In a vicinity of the SIM, slow and fast motions 



are locked in the image and null space of the thermodynamic projector P 



respectively. 



space of maps ^(0 [15] 



dm 

dt 



fm))-pfwt))- 



(ii) 



Rigorous proofs of existence and uniqueness of SIM, by the film equation ( llljl . were recently 
given for linear systems while the rationale behind the ( TIT]) is explained by means of 
a cartoon in the Fig. EK). Here, it is worth stressing that the above (llljl denotes a partial 
differential equation (PDE) whose unknown is a mapping (£) from a low dimensional 
reduced space H - £ G H - (also referred to as parameter space in the following) into the 
phase space U -tp EU . Therefore, readers should not get confused between stable stationary 
solutions of ( TIT]) (defining SIM as a mapping from H into {/) and single stationary states (or 
equilibrium states) of ([1]) (which satisfy the condition: f (ip ss ) = 0). 

For thermodynamically consistent systems ([1]) equipped with a potential G (thermody- 
namic Lyapunov function with respect to (JT])), MIM offers a projector whose construction is 
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based on the tangent space T and the gradient of the thermodynamic potential, dG/dip, at 
every point of SIM. This consistently imposes that the reduced dynamics P/ (?/>(£)) is dissi- 
pative. Explicit formulae for this thermodynamic projector are not necessary for the scope 



of this paper, and can be found in 



15] . Importantly, separation of motions in a vicinity of 



SIM is dictated by thermodynamic projector P, since it can be proved that slow motions 
along SIM are locked in the image, imP = T, whereas the null-space, kerP, spans the fibers 
of fast motions transversal to SIM (Fig. [Ib)) 16 1. 



Final 
of MIM 



y, a computationally advantageous realization is provided by a grid representation 



18] , where grid nodes in the phase space are defined by a discrete set of macroscopic 



variables, £, while finite difference operators are used to compute the tangent space at every 
node Thanks to locality of MIM constructions, we further make no distinction between 

manifolds and grids. 

i?emarfc-Consistent constructive methods of slow invariant manifolds rely upon efficient 
methods for solving the PDE ( TlTj) . As discussed below in section Hi I A 
finite difference schemes have been suggested in the literature 15 



towards this aim, 



20] (see also (EE 



Nevertheless, to the best of our knowledge, only explicit (or semi-implicit) schemes are 
available so far. Thus, due to hyperbolicity of the equation f JTTj) . its numerical solution is 
hindered by numerical instabilities (i.e. Courant type) 19], and no satisfactory solution to 
this issue has been suggested up to now. It is useful to stress that here we review the notion 
of film equation only for a better understanding of the present work. In fact, our suggestion 
toward the effective answer to the above problem is to avoid direct solution of ( TlTJ) (e.g. 
by finite difference schemes) in favor of its emulation, where the problematic term —Pf is 
not approximated with finite differences but mimicked by a redistribution step in terms of 
macroscopic variables (see section ITVl below) . 



A. Direct solution of the film equation 

A natural approach to the construction of SIM's is a direct numerical solution of the film 
equation f TTT]) starting with an initial (usually non invariant) manifold. For that, both the 
initial condition as well as implicit or semi-implicit schemes were developed. The simplest 
explicit scheme for solving the equation fTTT]) can be realized by iteratively refining each point 
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i/j of the initial manifold: ip + dip, 



# = r(/(^)-P/W), (12) 

with the time r being estimated according to the suggestions in jl8], where the scheme ([12]) 
is referred to as the relaxation method. It has been noticed [3] that the solution of the film 
equation of dynamics (TTTj) . similarly to hyperbolic partial differential equations for com- 
putational fluid dynamics (CFD) simulations, is hindered by severe numerical instabilities 
(see, e.g., the Courant instability 2l|). Furthermore, we notice that, unlike CFD, numerical 
solution of ( TTTj) comes with additional difficulties, due to an uncontrolled variation of the 
grid-node spacing. As a result, it is difficult to formulate an analog of the CFL (Courant 
- Friedrichs - Lewy) condition 2l| for (fl2|) . and the suppression of instability was only at- 



19] . In general, the latter 



tempted by an arbitrary decrease of the time r until convergence 
approach proves rather poor since the lack of convergence of ( f!2l) might not have numerical 
origin. In fact, there is no guarantee that the chosen number of reduced degrees of freedom 
q reveals sufficient in describing the asymptotic behavior of the dynamical system (CQ) in a 
given domain of the phase-space. For instance, in the case a higher number of reduced vari- 
ables are requested, the refinement of a q dimensional manifold by stable numerical schemes 
of (ITT]) is expected to fail anyway. The idea of adaptive dimension of SIM, formulated below 
in the section IIVBI is based on the latter observation. 

Finally, the construction of slow invariant manifolds by the solution of f fTT]) has been 
always attempted in the whole phase-space, by assigning a priori their dimension q somewhat 
arbitrarily. Such an approach, where the dimension q comes as external input into the 
problem, poses severe limitations to the accuracy of the reduced description and, most 
detrimentally, hinders the gaining of any better physical knowledge about it. Moreover, 
construction of high-dimensional invariant manifolds (q > 3) by the ( fTTT) is quite problematic 
and was never successfully accomplished up to now. 



IV. THE RELAXATION REDISTRIBUTION METHOD: RRM 

Toward the end of overcoming the above drawbacks, in this work, we introduce an ap- 
proach to model reduction, which allows for the construction of slow invariant manifolds 
with the dimension q adaptively varying from one region of the phase space to another. We 
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address thereby the fundamental issue of the minimal number of important (slow) variables 
which underlie the behavior of a complex dissipative phenomenon in a region of the phase 
space: A knowledge, emerging from the system and no longer imposed, is now gained. The 
latter is a challenging problem in Physics, and even in the classical cases, such as the re- 
duced description of the Boltzmann kinetic equation by a finite set of velocity moments of 
the distribution function (see, e.g., |22j), some essential questions remain open [111. l23l. \24\. 



Similarly, in chemical kinetics, several methods have been suggested |25M27l| for approxi- 
mating and parameterizing the SIM, however the choice of the minimal number of chemical 
coordinates (manifold parameters) is still debated. 

In the following, the key idea of our approach is to abandon an attempt of solving the 
film equation ffTTl) by numerical schemes such as ffT21 . in favor of a simulation of the physics 
behind this equation, in a spirit similar to Monte Carlo methods: As a consequence, a 
highly efficient construction of SIM with an embarrassingly simple implementation in any 
dimension is derived. 



A. Global formulation of RRM 

In order to introduce our method, we consider reaction kinetics and assume that a slow 
dynamics of (OQ) evolves on a g-dimensional SIM in the n- dimensional concentration space 
(this assumption will be relaxed in a sequel). Inspection of the right-hand side of fTTTT) 
reveals a composition of two motions: The first term, f(ip(£)), is the relaxation of the initial 
approximation to SIM due to the detailed kinetics, while the second term, —Pf(ip(^)) is the 
motion antiparallel to the slow dynamics. Let a time stepping St and a numerical scheme 
(e.g. Euler, Runge-Kutta, etc.) be chosen for solving the system of kinetic equations: All 
grid nodes relax towards the SIM under the full dynamics / during 5t. Fast component of / 
leads any grid node closer to the SIM while at the same time, the slow component causes a 
shift towards the steady state (see Fig. Eh.)). As a result, while keeping on relaxing, the grid 
shrinks towards the steady state (we term this a "shagreen effect" per de Balzac's famous 



novel 



2s 



- chagrin in French). Subtraction of the slow component therefore prevents the 
shagreen effect to occur, and it is precisely the difficulty in the numerical realization: explicit 
evaluation of the projector P on the approximate SIM does not always balance the effect of 
shrinking. This leads to instabilities, and results in a drastic decreasing of the time step. 
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FIG. 2: (Color online), a) The relaxation due to dU of a non-invariant manifold. Fast dynamics 
drives it toward the slow invariant manifold, whereas the concurrent action of the slow dynamics 
causes a shift toward the steady state (shagreen effect). On the contrary, relaxation due to the 
film equation (llip - (|12p allows movements only in the fast subspace. b) Relaxation Redistribution 
Method. The displacement in the slow subspace, generated during relaxation, is annihilated by a 
redistribution step in the parameter space. 

The key idea here is to neutralize the slow component of motion by a redistribution 
of the points on the manifold after the relaxation step (see Fig. Wp)). For the sake of 
presentation, we assume that macroscopic parameters are given by a set of q linear functions 
b = {bi, . . . , b q } such that b\ (if)) = . . . , b q (ip) = £ 9 . Let £ = (£ 1 , . . . , £ 9 ) be a generic node 
of a fixed grid S in the parameter space S, and the g-dimensional slow invariant grid (SIG) 
in the phase space U is initialized: Q m = if) m (^) (that is, the initial SIG is the collection of 
nodes if) m = if) m (£,), After the relaxation step, all the nodes ip m have moved to new 

locations, ip m — > and we denote £ R = b(ip R ) the values of the macroscopic parameters 
corresponding to the relaxed nodes i/j r . 

It is worth stressing that by parameter space here we mean the low dimensional macro- 
scopic space S whose dimension is q << n. Hence, an arbitrary grid S is defined by a 
mapping, if; (£), on a subspace of H into the phase-space U (of dimension n). 
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For example, the forward Euler scheme used below gives 

^=ft)+«/rt). (13) 

With this, also the nodes of the grid S shift by an amount <5£ = b(^ R ) — £ due to the slow 
component of motion. The redistribution of the nodes ip R back to the fixed grid S simulates 
the subtraction of the slow motion from the relaxation step, and is done as follows: For 
each £ G S, we consider a g-simplex S q (in U) with q + 1 vertices ^> R , , . . . , ^ R such that 
£ is inside the macroscopic projection of S q , the simplex H q (in S) formed by the vertices 
£ R = b(ipf), £ R = b(ipf-), . . . , £ R = b(ip R ). The updated (relaxed-and-redistributed) grid 
qRR j g constructed by a linear interpolation of the vertices of the simplex S q : 

^ RR = ( 1 - E w ) ^ + E ( 14 ) 

V i=i / i=i 

where the weights lOj are so chosen as to satisfy the redistribution condition, 

^# RR ) = t- (15) 
This amounts to solving a q x q linear system, 

iy>jW?) - bM)]*>j =C~ bM)- (16) 
i=i 

The above procedure is supplemented by the boundary conditions applied at the edges of the 
grid: Grid nodes at the boundary ip^ are reconstructed by extrapolation after the relaxation 
step. Formula (III"]) is used where ip RR = ipb £ S q is located in the vicinity of a simplex S q 
with vertices ip R , ijj R , . . . , ip R . In general, S q can be chosen in such a way that its vertices 
are the relaxed states of the initial nodes ip™, ip™, . . . , ip q n with ip^, = ip™. 

Thus, after the redistribution step, the initial grid is refined towards the invariant grid. 
The procedure is then iterated, whereas each relaxation step is altered by the redistribution 
step, in which the slow motion is subtracted by stretching the macroscopic variables to the 
nodes of the initial grid S. 

We notice that, on SIM, movements due to the vector field / occur along the manifold 
itself, thus the effect of the relaxation is entirely counterbalanced by the subsequent redis- 
tribution on the SIM. It is worth stressing that this observation holds for every invariant 
manifold (not necessarily SIM). Nevertheless, numerical evidences clearly show that an arbi- 
trary invariant manifold Qi nv is an unstable solution of the above dynamics, and refinements 
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starting from Vti nv converge toward the SIM, which instead turns out to be a stable solution. 
As a result, slow invariant grids are stable stationary solutions of the described procedure, 
here termed relaxation redistribution method (RRM). Once the invariant grid is constructed, 
the reduced dynamics for variables £ is defined as 

§ = W RR (£))- (17) 

In other words, the suggested RRM enables to provide the reduced system (|17[) . written in 
terms of a significantly smaller set of variables £, with a closure. 

Note that, upon the global construction of SIG, computations deliver a discrete set of 
linked states ip RR {!;), i n a vicinity of the corresponding slow invariant manifold. Here, grid 
nodes are termed "interconnected" because we assume that for any arbitrary node it is pos- 
sible to identify all its nearest neighbors. Moreover, interconnectivity enables one to easily 
proceed with analytical continuation of the above slow invariant grid, and thus to the calcu- 
lation of the right-hand side of ( |T7|) for any set of variables £. To this end, for simplicity, here 
we adopt multi-linear interpolation, which posses the advantage to automatically fulfill the 
linear conservation constraints ([8]). For further details on multi-linear interpolation of grids, 



the interested reader is delegated to 29j] . On the other hand, if the local construction of SIG 
is implemented, a closure for (1171) is computed when needed and no analytical continuation 
of the grid is requested. In the latter case, in order to speed up the computations, smart 
methodologies for data storage and retrieval can be used and are readily available from the 
literature (see, e.g., the ISAT method in jsoj). 

Finally, note that while the redistribution step seems " natural" from the numerical stand- 
point of discretizing the above macroscopic equation (fT7l) on a fixed grid S, the feature 
recognized here is that it is precisely the subtraction of the slow component of the motion 
in the film equation (II ip . which circumvents the question about explicit evaluation of slow 
motions in the course of the SIM construction. 

In order to test the RRM, we first consider a simple benchmark suggested by Davis and 
Skodje (DS) jsi| (a two-dimensional system with a one-dimensional SIM known in a closed 



analytical form). The DS system 3JJ| consists of two equations 



dx/dt = f x (x) = -x, ^ 
dy/dt = f y (x, y) = -jy + [(7 - 1) x + -fx 2 ]/(l + xf , 7 > 
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FIG. 3: (Color online). The Davis-Skodje system 31[. Two different initial grids are refined 
using the forward Euler scheme for the relaxation (St = 10 -2 ). Results after 50 RRM iterations 
are reported (refined grids) with 7 = 50. Triangles show an intermediate step (after two RRM 
iterations) starting from the initial smooth grid. 



has unique stable steady state x = y = 0, and a one-dimensional SIM, y = x/(l + x). Here, 
when 7 >> 1, due to a significant separation between of time scales of the two variables x and 
y, all solution trajectories of DS system exponentially decay to the SIM (see Ref. [31|). In the 
above notation, ip = (x, y) T , and we define the slow variable as £ = x, that is, b = (1, 0) and 
b(ip) = (1, 0) ■ (x, y) T = x. The RRM is initialized with the grid represented by the collection 
of points {(x r , y m (x r ))}, where x r are distributed evenly in the interval x G [0, x\>]- Upon the 
relaxation step, the grid points are shifted to new locations {(x r , y m (x r ))} —> {(x R , y R (x r ))} 
with x R = (1 — St)x r , y R = y m (x r ) + 5tf y (x r ,y m (x r )). Choosing the interval Si = [x^, xf] 
for each point x r such that x r G £1, the redistribution ( fT4l) gives 



V 



RR 



[X 



%o)yf 



xf X R 



(19) 



while x RR — x r , by the condition (|T5l) . For the boundary node at x^ we set y R = y R and 
for y R = y R (xb_i), with x^-i G S being the nearest neighbor of x^. In Fig. [31 the local grid 
step at the boundary is: 5xb = Xb — Xb-i = 6 — 5.88. 
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The RRM was performed for a variety of initial grids (initialized with different functions 
y m (x)), with different spacing and with various choices of the simplex. Independent of these 
variations, the RRM iterations converged stably to the analytical SIM of the DS system. 
Results are presented in Fig. [3] for two different initial grids, a regular linear (y m = ax, 
a = 0.25) and a randomly generated grid (for each value x r G [0, 6] a random number 
y m (x r ) e]0, 1[ is assigned by a linear congruent generator) with the intervals Si chosen as 
X\—Xo = 0.12. Convergence to SIM is even striking in Fig. [3jgiven the fact that both initial 
grids are far from SIM. 

Thus, convergence of the RRM iterations confirms the existence of a reduced description 
with a fixed number of degrees of freedom q (existence of g-dimensional slow invariant 
manifold) . On the contrary, no convergence in RRM indicates that more degrees of freedom 
are needed to recover the detailed system dynamics. This concept shall be used below for 
adaptively choosing the invariant grid dimension. 

Both construction and usage of a global reduced description soon become impracticable 
as the dimension q increases. In fact, computing and storage of high dimensional SIM's may 
be problematic already at q > 3. Above all that, data retrieval by interpolation on such 
large arrays is computationally intensive, and sometimes full construction of manifolds can 
be useless: For example, in combustion applications, regions with high a concentration of 
radicals are unlikely to be visited. 

Remark-In general, when using model reduction techniques, such as the RRM method, 
slow and fast subspaces are not known in advance. In fact, this kind of information is what 
we get at the end of the process. Invariant grids constructed by the suggested RRM are 
finally located in the slow subspace (regardless of the choice on the parameterization). The 
fast subspace can be thereafter reconstructed by adopting e.g. the notion of thermodynamic 



projector (see, e.g., Refs. 15|, [l6|, ll9|). On the other side, concerning the parameterization 



choice, we notice that (as stressed in the conclusions I VI I) there are no universal recipes, 
and it specifically depends on the physical phenomenon we are dealing with. In general 
good macroscopic variables can be found in the literature: For instance, in the case of 
the Boltzmann equation typical macroscopic parameters are the velocity moments of the 
distribution function, whereas for chemical kinetics we can use spectral variables as done for 
the example in section V. Alternatively, in the latter case, typical slow variables can also be 
adopted (see, e.g., the RCCE parameterization in |25|). 



15 



B. Local formulation of RRM 
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FIG. 4: (Color online), a) Relaxation redistribution method: Local formulation. Only a small 
patch of the SIM is constructed. After refinement, the coordinates of the pivot provide the reduced 
system (|17j) with a closure, b) Simplexes can be conveniently adopted for a patch-wise description 
of the SIM in any dimension. 



Importantly, the RRM allows for a straightforward local formulation, where only small 
patches of the slow invariant grid are initialized and refined. Let £ = (f 1 , . . . , £ 9 ) and 
the procedure is initialized with a simplex S q where the pivot ip m = ip m {!;) is linked to 
q secondary nodes i/i™ = V' m (£i); ip 1 ? = ?/>(6j) in a neighborhood of ip such that £j = 
. . . , £ l + 5£\ . . . , £ 9 ) , with 6£ l being a small deviation of the z-th macroscopic variable. 
A sequence of relaxation and redistribution steps is applied to the vertices of S q in any 
dimension q: This realizes indeed the simplest instance of the RRM, 



RR 



(20) 



i=i 



i=l 



while the weights Wi are found from the redistribution (anti-shagreen) condition (fl~5l) : 
b(vjj RR ) = £. Refinements end as soon as a norm of the total displacement of the pivot 



at the nth RRM iteration, 



RR 

(n+1) 



, becomes sufficiently small compared to 



the displacement caused by the relaxation alone, 

1(3 



RR 



(n+1) 



>) 



Setting an upper limit to both the number of refinements N and the tolerance e such 
that: 

W)\ / \S$$\ < e, (21) 

the local RRM can be adaptively performed starting with q = 1. If the latter requirements 
are not fulfilled, the dimension is updated to q = 2 and the procedure repeated. Upon 
convergence with some q = q, a closure of the reduced system (JTTj) is provided by the coor- 
dinates of the pivot. It is worth stressing that the above convergence criterion (}2~Tj) is based 
on the value assigned to the tolerance e and number of refinements N. However, the latter 
quantities can be properly set upon an independence study with respect to the manifold 
dimension. Namely, in the same spirit of grid independence studies of fluid dynamics simu- 
lation results, the independence of the manifold dimension q on e and N can be verified by 
repeating the calculations with smaller tolerances and larger number of refinements. In this 
sense, the local RRM fully alleviates any assumption about the dimensionality of SIM, the 
local dimension is found automatically and if no reduced description is possible at all, no 
convergence at any q < n will clearly indicate this. 

Finally, for systems supported by a Lyapunov functions G (such as the kinetic equations 
(jSJ)), a convenient (but not the only possible) initialization of the RRM procedure (construc- 
tion of the initial pivot and secondary nodes) for dissipative systems can be accomplished 
by means of the notion of quasi equilibrium manifold (QEM). In this respect, an approxima- 
tion of the g-dimensional SIM can be obtained by minimizing the function G under q linear 
constraints in addition to the element conservation laws (jHJ): 

G — > min 

MV) = e, i = l,...,<Z (22) 
Cijj = const. 

where, in the case of chemical kinetics, the function G is a thermodynamic potential (i.e., 
entropy, Gibbs free energy, etc.) as discussed in the section Hi Al It is worth stressing that 
the idea of using extrema of potentials, for providing reduced description with a closure, 



dates back to the work of Gibbs 



321 ] . From then on, this notion has been adopted is several 



areas such as the kinetic theory of gases [15j, |33[, or detailed combustion mechanisms [25 1. 



However, we stress that the latter approximations often provide with a poor description of 



the corresponding SIM 27|, |34j, thus they are used here only for initializing the RRM. 
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Below, following the suggestion in 35[], we make use of spectral variables £ l = bi(ip) 
obtained by the inner product between the state i/j and the parameterization vectors bi, 
which are the left eigenvectors of the Jacobian J = df /dip at the steady state, corresponding 
to non-zero eigenvalues A« and numbered in the order of increase of |Aj|. The latter is referred 
to as spectral quasi equilibrium parameterization. The pivot ip* = (ipi, ...,ip n ) of the initial 
simplex S q is defined as the quasi-equilibrium point |l5|], corresponding to £ = (£ 1 , . . . , 
and calculated by solving the problem fl22l) . To this end, the (1221) is equivalent to the global 
minimization problem of a Lagrange function G: 

G = G + J2i b *W- C] A, + ACty, (23) 

i=l 

with Aj, A being a set of Lagrange multipliers. We notice that, efficient tools for the solution 
of fl22|) are also available (see, e.g., STAN JAN 36|). Secondary nodes ip k of the simplex 



inear expansion of the minimization problem about the 



3| : ^ = ^* + YJiZfKpi, with (px, . . . p n „ d ) and 5 k 



can be conveniently calculated by 
quasi-equilibrium as suggested in 
(o~l, ^k~ d ) being a vector basis spanning the null space of C and the solution of a linear 
algebraic system 

' Y.Zt (tjH*Pi) = -VG%-, j = 1, . . . , n - d - q, 

YllZi (hpi) Si = o, 



(24) 



k TJlJ {bqPi) $i = 0. 

The vector basis (t\, . . ,t n _ g _d) spans the kernel of the linear space defined by the vectors 
bi and the rows of the matrix C in (jSJ), H* = [d 2 G / dipidipj] and VG* = [dG/dipi] are the 
second derivative matrix and the gradient of the function G at the pivot respectively, while 
£fc defines the length of the edge of the simplex S 9 along the k-th direction. 



V. ILLUSTRATION: DETAILED HYDROGEN-AIR MIXTURE 

Here, we considered the autoignition of hydrogen- air mixture at stoichiometric proportion, 
reacting according to the realistic detailed mechanism of Li et al. j^j, where nine chemi- 
cal species and three elements participate in a complex reaction dictated by twenty-one 
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FIG. 5: (Color online). Heterogeneous slow invariant manifold of hydrogen-air combustion mecha- 
nism by local RRM. Three-dimensional projection of the six-dimensional phase space onto spectral 
variables (see text). The two-dimensional patch ("kite", triangles) is tight by a one dimensional 
" thread" (line) to the zero-dimensional equilibrium and merges with the three-dimensional " cloud" 
(tetrahedra) . Legend: mass fraction of OH. Explicit Euler scheme with 5t = 5 x 10 -8 [s] was used 
for the relaxation of simplex. RRM convergence criteria: N = 2000, e = 10 -4 . 

reversible elementary steps (J2J) (this mechanism is universally used in turbulent combustion 
simulations {4], and details for this case are discussed in the Appendix |A]) . Time evolution 
of species concentration is governed by ([3]), and it is supplemented by the condition for the 
reactor temperature, which stipulates the conservation of the mixture enthalpy (adiabatic 
reactor) : 

9 

h = Y^ h i y i = 1000[JfeJ/Atf]. (25) 

i=l 
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Furthermore, the pressure of the mixture is fixed (p = l[atm]), and the mass fraction (Yi) 
of an arbitrary chemical species i can be expressed in terms of the corresponding molar 
concentration (q) by means of the following relationship: 

Y t = " lUt ■ (26) 

Fig. [5] shows a projection of the heterogeneous SIM (i.e. with a varying dimension in 
the phase-space), onto the subspace £ 1 , £ 2 , £ 3 , Yqh, constructed by the local RRM where 
one-, two- and three-dimensional patches are clearly visible. Here, the variables £ 1 , £ 2 , £ 3 are 
chosen according to the spectral quasi equilibrium parameterization, where = bi(ip) = 
with hi denoting the three slowest eigenvectors of the Jacobian matrix J = df /dip at the 
steady state, whereas the RRM is initialized as discussed above in the text with the potential 
G computed on the basis of the mixture-averaged entropy fl9]). Interested reader may find 
fill details on the computation of G and its derivatives (VG and H* requested in (I2~3| ) in 
37|. Results in terms of basic variables (i.e. concentrations of species) can be obtained upon 
a post-processing of the spectral variables which amounts to a linear transformation. 

A typical problem, where dynamics evolves along a cascade of slow invariant manifolds 
with progressively lower dimensions, is the auto-ignition of a fuel-air mixture. In Fig. EJ 
solution of the reduced system (fT7|) . supplemented with a closure by the local formulation of 
RRM, is compared with the integration of the detailed reaction mechanism. Results are in 
excellent agreement for all the chemical species and the temperature. Note that, although 
one- and two-dimensional SIM's are able to recover the most of the dynamics of major 
species and of the temperature, the minority species (such as radicals HO2 and H 2 02) do 
require high dimensional manifolds (q > 3) to be correctly predicted. 

For the sake of clarity, we outline below the steps leading to the computation of a q- 
dimensional closure corresponding to a macroscopic state £ = £ 9 ), by the local RRM 

for the above kinetic system: 

1. Set up the initial SIM dimension (e.g., q = 1); 

2. Set up a convergence criterion (12T1) . and the maximal number of iterations iV; 



3. Compute the initial coordinates of the pivot ip*, which amounts to solving a non linear 
algebraic system, VG = 0, e.g. by Newton-Raphson iterations; 
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4. Compute q secondary nodes ip k by the linear algebraic system ( 124)) ; 

5. Update the coordinates of both the pivot and secondary nodes by the RRM equation 
(EQJ); 

6. Check convergence; 

7. if no convergence is achieved after N iterations, then update the SIM dimension q = 
q + 1 and go to 3; 

8. exit. 

The above illustration demonstrates that the suggested RRM method is able to accu- 
rately recover the dynamics of a complex system. Moreover, here we adopted the automatic 
criterion (I2TT) to choose the number of reduced degrees of freedom (macroscopic important 
variables), which are strictly needed to reproduce the phenomenon under study. The latter 
features make the RRM, on one side, a pretty useful tool for the efficient computation of 
large dissipative systems. Most importantly, on the other side, it enables to gain a better 
physical understanding about a complex phenomena by addressing the issue of its minimal 
description. 



VI. CONCLUSION 



To conclude, we addressed here the fundamental problem of the minimal description of a 
complex dissipative system which is a challenging issue in Physics. Our approach is based 
on a simulation (instead of a solution) of the fundamental film equation of dynamics (TTTT) . 
We stress that it is the RRM realization which is able to unfold the full power of the Method 
of Invariant Manifold which was not possible before, such as the adaptive construction of 
high dimensional manifolds (i.e. q > 3, with q varying from a region of the phase-space to 
another). On the practical side, RRM is pretty simple as it is based on a direct integration of 
the film equation plus redistribution. The key point realized in this paper is that the latter 
simulates subtraction of slow motion from the film dynamics, the step which is hard to 



control in more conventional approaches to the film equation 



ljj. In that respect, the RRM 



is similar in its spirit (but certainly not in the implementation ) to other successful simulation 



strategies such as the Direct Simulation Monte Carlo method 
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which replaces the solution 
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FIG. 6: (Color online). Auto-ignition of homogeneous stoichiometric mixtures of hydrogen and 
air: Histories of the temperature and of the mass fraction of chemical species. Line: detailed 
reaction; Symbol: local RRM method by adaptively following a cascade of reduced models of 
various dimension q: q = 5 (square), q = 4 (cross), q = 3 (diamond), q = 2 (star), q = 1 (circle) 
and q = (steady state). 

of the Boltzmann equation by a stochastic simulation of "collisions" . We stress that, suitable 
macroscopic variables depend on the specific phenomenon (e.g. velocity moments of the 



distribution function for describing gas kinetics [22|, |39(). The methodology developed in 
this paper addresses the general problem of minimal macroscopic description by letting the 
system decide how many important variables are to be considered. 

Examples presented above convincingly show that RRM achieves all the objectives set for 
obtaining the accurate reduced description, whereas the resulting adaptively reduced models 
reveal new physical knowledge of a complex dissipative system (i.e. its minimal description), 
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and can be used for its computationally efficient simulation. We should stress that fully 
adaptive construction of heterogeneous slow invariant manifolds as in the case of hydrogen- 
air mixture is difficult if at all possible with any other model reduction technique 7| • Finally, 
while we focused on the important class of dissipative systems arising in combustion, we look 
forward to generalization of the above technique of simplification to other dissipative systems 
such as master and Fokker-Planck equations and other complex dynamics. 
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Appendix A: Detailed reaction for hydrogen and air 



In the table [H we report the list of all reaction steps involved in the combustion mechanism 
for hydrogen and air adopted in section [V] where n = 9 species (H 2 , N 2 , H, O, OH, 2 , 
H2O, HO2, H2O2) and d = 3 elements (H, O, N) are involved in r = 21 elementary 
reversible steps. The system of kinetic equations is formulated according to the ([3]) and (j4]), 
where the reaction constant k^ of the s-th step is determined by the Arrhenius law ([5]) with 
the coefficients A s , n s and Ea s from table [H In the following, the symbol M represents an 
additional species, whose concentration cm denotes a weighted sum of the concentration of 
all species (third-body reaction): 

n 

cm = y^jCu ( A1 ) 
i=i 

cjj being the third-body efficiencies. In the reactions N. 5, 6, 7 ,8, it is adopted an 2 o — H-0, 
aH 2 = 1-5, and a« = 1 for all other species. Finally, the steps N. 9 and 16 are typical fall- 
off' reactions, where the reaction constant remarkably depends on the mixture pressure. 
In this case, fc+ and k$ are the reaction constants in the high- and low-pressure limit, 
respectively, and the reaction constant reads: 

fc+ = k+ FP r f (1 + P r ) , (A2) 



with P r = k^CM/k^, and F given by the Troe function (see 40| for the details). In particular, 
in the reaction step N. 9 the third-body efficiency are clh 2 o — 10, ao 2 = —0.22, in the reaction 
step N. 16 an 2 o = 11, < 1 h 2 =1.5, whereas in both cases a« = 1 for the rest of the species. 
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TABLE I: Detailed 772-air reaction mechanism. Units are cm 3 , mol, sec, Kcal and K. a Troe 
parameter is: 0.8. 6 Troe parameter is: 0.5. 
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